function [LQi, vHat, Lrij_out,Lrji] = softldpc1(Lci, H,nozinrow,nozincol, iteration) %#codegen
% Log-domain sum product algorithm LDPC decoder
%
%  rx        : Received signal vector (column vector)
%  H         : LDPC matrix
%  N0        : Noise variance
%  iteration : Number of iteration
%
%  vHat      : Decoded vector (0/1)
%
%
% Copyright Bagawan S. Nugroho, 2007
% http://bsnugroho.googlepages.com
NMS = 1;
[M, N] = size(H);

% Prior log-likelihood. Minus sign is used for 0/1 to -1/1 mapping
vHat = zeros(1,N);
LQi=Lci;
limit_Lr = 0.8;
% Initialization
Lrji = zeros(M, N);

alphaij = zeros(M,N);
betaij = zeros(M,N);
Lrij_out = zeros(M,N);
% Asscociate the L(ci) matrix with non-zero elements of H
Lqij = H.*repmat(Lci, M, 1);%%%%%%换成下面的%+1e-80这种做法可能有问题
%Lqij = Lqij_in;

% Iteration
for n = 1:iteration
    
    %fprintf('Iteration : %d\n', n);
    
    % Get the sign and magnitude of L(qij)

     if(n == 1)
      alphaij = sign(Lqij + 1e-80);
      betaij = NMS.*abs(Lqij);
     end
    
    % ----- Horizontal step -----
    for i = 1:M

        c1 = nozinrow{i};
 %        tempArray = betaij(i,c1);
        [temp_sortted,~] = sort(betaij(i,c1));
        prodofrow = prod(alphaij(i, c1));
        
%2021.11.9
%        prodoftanh_Lqji = prod(tanh_Lqji(i,c1));
%2021.11.9
        

        % Get the summation of Pi(betaij)
        for k = 1:length(c1)
            if( temp_sortted(1) == betaij(i,c1(k)))
                sumOfPibetaij = temp_sortted(2);
            else
                sumOfPibetaij = temp_sortted(1);
            end
            
%             sumOfPibetaij = 0;
%             prodOfalphaij = 1;
%             tempArray = betaij(i,c1);
%             tempArray(k) = [];
            
            % Summation of Pi(betaij)\c1(k)
            % sumOfPibetaij = sum(Pibetaij(i, c1)) - Pibetaij(i, c1(k));
%             sumOfPibetaij = min(tempArray);
%             prodOfalphaij = prod(alphaij(i, c1))*alphaij(i, c1(k));
              prodOfalphaij = prodofrow * alphaij(i, c1(k));
            
            % Update L(rji)
              Lrji(i, c1(k)) = prodOfalphaij*sumOfPibetaij;
%2021.11.9
%        divide_prodoftanh_Lqji = prodoftanh_Lqji/tanh_Lqji(i,c1(k));
%        Lrji(i,c1(k)) = atanh(divide_prodoftanh_Lqji);
%        if((Lrji(i,c1(k)))>limit_Lr)
%            Lrji(i,c1(k)) = limit_Lr;
%        else if(Lrji(i,c1(k))<-1*limit_Lr)
%                Lrji(i,c1(k)) = -1*limit_Lr;
%            end
%        end
%        Lrji(i,c1(k)) = 2*Lrji(i,c1(k));
%2021.11.9
           
        end % for k

    end % for i
    
    % ------ Vertical step ------
    for j = 1:N
         %test
         if(j== 1216) 
             a=1;
         end
        % Find non-zero in the row
%         r1 = find(H(:, j));
        r1 = nozincol{j};
        Lrjisum = sum(Lrji(r1, j));
        
        for k = 1:length(r1)
            
            % Update L(qij) by summation of L(rij)\r1(k)
%             Lqij(r1(k), j) = Lci(j) + sum(Lrji(r1, j)) - Lrji(r1(k), j);
              Lqij(r1(k), j) = Lci(j) + Lrjisum - Lrji(r1(k), j); %+1e-80这种做法可能有问题
              Lrij_out(r1(k), j) = Lrjisum - Lrji(r1(k), j);
%               if(Lqij(r1(k),j) == 0)%%%%%%%%%%%%%%%%%%
%                   alphaij(r1(k),j) = 1;%%%%%%%%%%%%%%%%%%%%%
%               else%%%%%%%%%%%%%%%%%%%%%%%%
            alphaij(r1(k),j) = sign(Lqij(r1(k),j)+ 1e-80);
%               end%%%%%%%%%%%%%%%%%%%%%%
            betaij(r1(k),j) = NMS.*abs(Lqij(r1(k),j));
%2021.11.9
% %tanh(0.5lqji)
%         lqji(r1(k),j) = 0.5.* Lqij(r1(k),j);
%         tanh_Lqji(r1(k),j) = tanh(lqji(r1(k),j));
%2021.11.9
        end % for k
        
        % Get L(Qi)
%         LQi(j) = Lci(j) + sum(Lrji(r1, j));
        LQi(j) = Lci(j) + Lrjisum;
       % fprintf('LQi: %d\n',LQi);
        % Decode L(Qi)
        
        
    end % for j
    for fin=1:N
    if LQi(fin) < 0
                vHat(fin) = 1;
             else
                vHat(fin) = 0;
    end
    end
%     if(H*vHat'==0)
%         break;
%     end  

end % for n